/*
 * Geotoolkit.org - An Open Source Java GIS Toolkit
 * http://www.geotoolkit.org
 *
 * (C) 2018, Geomatys.
 *
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation;
 * version 2.1 of the License.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Lesser General Public License for more details.
 */
package org.geotoolkit.util.grid;

import java.util.Arrays;
import java.util.Spliterator;
import java.util.function.Consumer;
import java.util.function.DoubleUnaryOperator;
import static org.geotoolkit.util.grid.GridTraversal.isNearZero;
import static org.geotoolkit.util.grid.GridTraversal.EPSILON;
import static org.geotoolkit.util.grid.GridTraversal.toVector;

/**
 * Find intersection points between a segment and a regular grid whose cells sides are of length 1. Works in
 * N dimensions.
 *
 * @implNote To find intersection with the grid borders, we do not use equation of intersection between the vector
 * of the segment and the plans of the gride edges. Instead, we use colinearity properties. In fact, at each step, we
 * search a vector overlapping the segment, but stopping at the next encountered grid edge. That means that for a
 * segment AB, we search the shortest vector AC, where:
 * <ul>
 * <li>AB = kAC, where k is a positive (same direction), non null value.</li>
 * <li>C is the closest distinct point to the C-1 point (which is C point computed at previous step), with at least one
 * of its ordinate is an integer.</li>
 * </ul>
 *
 * Optimisations: There's a lot of chances that a more efficient algorithm exist. For now, the complexity is O(2*d*n),
 * where d is the number of dimensions of input points, and n is the number of intersections between the segment and
 * the grid.
 *
 * @author Alexis Manin (Geomatys)
 */
public class MultiDimensionMove implements Spliterator<double[]> {

    /**
     * Start of the segment to browse along the grid
     */
    private final double[] startPoint;

    /**
     * Last point of the segment.
     */
    private final double[] endPoint;

    /**
     * Vector representing the segment to browse.
     */
    private final double[] segmentVector;

    /**
     * The latest point generated by this spliterator.
     */
    private final double[] previousPoint;

    /**
     * For each axis, it embed the operator in charge of moving the corresponding ordinate.
     */
    private final DoubleUnaryOperator[] moveDirs;

    /**
     * Create a new spliterator to find intersections between a segment AB and a regular grid.
     *
     * @param startPoint The start point of the segment to analyze (A).
     * @param endPoint The end point of the segment to analyze (B).
     */
    public MultiDimensionMove(double[] startPoint, double[] endPoint) {
        if (startPoint.length != endPoint.length) {
            throw new IllegalArgumentException(String.format(
                    "Dimension of given points does not match: start is %dD, but end is %dD",
                    startPoint.length, endPoint.length
            ));
        }

        this.startPoint = Arrays.copyOf(startPoint, startPoint.length);
        this.previousPoint = Arrays.copyOf(startPoint, startPoint.length);
        this.endPoint = Arrays.copyOf(endPoint, endPoint.length);

        // create the vector representation of our segment. In the meantime, we also prepare the operators in charge of
        // moving ordinate values.
        moveDirs = new DoubleUnaryOperator[startPoint.length];
        segmentVector = new double[startPoint.length];
        for (int i = 0; i < segmentVector.length; i++) {
            segmentVector[i] = endPoint[i] - startPoint[i];
            if (isNearZero(segmentVector[i])) {
                moveDirs[i] = null; // Ignored dimension. No move in this direction.
            } else if (segmentVector[i] < 0) {
                moveDirs[i] = GridTraversal::floorOrDecrement;
            } else {
                moveDirs[i] = GridTraversal::ceilOrIncrement;
            }
        }
    }

    @Override
    public boolean tryAdvance(Consumer<? super double[]> action) {
        double minDist = Double.POSITIVE_INFINITY;
        int idxOfMinDist = -1;
        for (int i = 0; i < moveDirs.length; i++) {
            DoubleUnaryOperator mover = moveDirs[i];
            if (mover != null) {
                double move = mover.applyAsDouble(previousPoint[i]);
                if (segmentVector[i] >= 0 ? move >= endPoint[i] : move <= endPoint[i]) {
                    move = endPoint[i];
                }

                move -= previousPoint[i];

                if (!isNearZero(move)) {
                    // find if the overall norm of the resulting vector is shorter than the previous one
                    final double moveRatio = move / (segmentVector[i]);
                    assert moveRatio >= 0.0;
                    if (moveRatio < minDist) {
                        minDist = moveRatio;
                        idxOfMinDist = i;
                    }
                }
            }
        }

        if (idxOfMinDist < 0) {
            return false;// We're already on the end point
        }

        for (int i = 0; i < previousPoint.length; i++) {
            previousPoint[i] += minDist * segmentVector[i];
        }

        assert areColinear(segmentVector, toVector(startPoint, previousPoint));

        action.accept(Arrays.copyOf(previousPoint, previousPoint.length));
        return true;
    }

    @Override
    public Spliterator<double[]> trySplit() {
        // TODO: it is splittable, but require non-trivial analysis: we must perform the same kind of analysis as in
        // tryAdvance, but for approximately half the overall distance of input segment.
        return null;
    }

    @Override
    public long estimateSize() {
        // TODO : how to (roughly) estimate the number of crossed grid points ? Segment norm is not efficient enough,
        // due to the distorsion of distances over big dimensions.
        return Long.MAX_VALUE;
    }

    @Override
    public int characteristics() {
        return ORDERED | DISTINCT | NONNULL;
    }

    static boolean areColinear(final double[] u, final double[] v) {
        double[] coefs = new double[u.length];
        for (int i = 0 ; i < u.length ; i++) {
            double uv = v[i] - u[i];
            if (isNearZero(uv)) {
                coefs[i] = Double.NaN;
            } else if (isNearZero(u[i])) {
                return false;
            } else {
                coefs[i] = v[i] / u[i];
            }
        }

        for (int i = 1 ; i < coefs.length ; i++) {
            if (coefs[i - 1] != coefs[i] && (coefs[i-1] - EPSILON > coefs[i] || coefs[i-1] + EPSILON < coefs[i])) {
                return false;
            }
        }

        return true;
    }
}
